Compensating for space and slowness/angle blurring of reflectivity

ABSTRACT

Systems and methods for compensating for spatial and slowness or angle blurring of plane-wave reflection coefficients in imaging. A wave field may be determined at a reference depth proximate to a reflector for a shot record. A receiver-side blurring function may be determined at the reference depth. An aggregate blurring function may be constructed based at least partially on the source wave field and the receiver-side blurring function. A plane-wave reflection coefficients may be determined based at least partially on the aggregate blurring function.

BACKGROUND

Seismic imaging may be used to delineate geological interfaces in asubsurface domain. Seismic imaging may begin with the acquisition ofcontrolled-source reflection data (e.g., marine data). In areas wherethe subsurface is laterally homogeneous, the type of seismic imagingused may be referred to as “time imaging.” In areas where the subsurfacehas greater complexity, the type of seismic imaging used may be referredto as “depth imaging.”

In addition to delineating geological interfaces in the subsurface, thematerial property contrasts at the interfaces may be useful, inparticular for interfaces at suspected oil and/or gas reservoirs.Amplitude versus offset or angle (“AVOA”) may be used to estimate theseproperty contrasts (e.g., in combination with time migration).

Applying AVOA may be a challenge for more complicated environments wheredepth imaging is used. For example, imaging below highly-refractive saltor basalt may be prone to illumination problems that may be quantifiedprior to AVOA inversion.

What is needed is an improved system, method, and resulting workflowsfor compensating for space and slowness/angle blurring of reflectivity.

SUMMARY

Embodiments of the disclosure may provide a method for compensating forspatial and slowness/angle blurring of plane-wave reflectioncoefficients in imaging. The method may include determining a sourcewave field at a reference depth proximate to a reflector for a shotrecord. A receiver-side blurring function may be determined at thereference depth. An aggregate blurring function may be determined basedat least partially on the source wave field and the receiver-sideblurring function. A plane-wave reflection coefficient may be determinedbased at least partially on the aggregate blurring function.

In an embodiment, the plane-wave reflection coefficient may further bebased at least partially on an interaction of the aggregate blurringfunction with a reflection operator containing the plane-wave reflectioncoefficient.

In an embodiment, the plane-wave reflection coefficient may bedetermined based at least partially on the aggregate blurring functionin a space-time domain at a fixed depth proximate to the referencedepth.

In an embodiment, the plane-wave reflection coefficient may bedetermined in a domain of lateral position and depth at a fixed time.

In an embodiment, a calibrated extended image gather may be constructedproximate to the reference depth based at least partially on migrationdata and a calibration field including a band limit of the receiver-sideblurring function.

In an embodiment, the reflector may be a dipping reflector with respectto a coordinate frame of the source wave field and the calibratedextended image gather.

In an embodiment, the calibrated extended image gather and the aggregateblurring function may be transformed to a space-frequency domain at afixed depth. A reflection operator may be obtained in thespace-frequency domain by matrix inversion.

In an embodiment, the reflection operator may be transformed from thespace frequency domain into the plane-wave reflection coefficient by aFourier transform.

In an embodiment, the aggregate blurring function may be transformedinto a smearing function in a slowness or angle domain of the plane-wavereflection coefficient.

In an embodiment, the smearing function in the slowness or angle domainmay be determined in a space-time domain at a fixed depth proximate tothe reference depth.

In an embodiment, the smearing function in the slowness or angle domainmay be determined in a domain of lateral position and depth at a fixedtime.

In an embodiment, the calibrated extended image gather may betransformed into the slowness or angle domain. The plane-wave reflectioncoefficient may be determined based at least partially on the smearingfunction in the slowness or angle domain.

In an embodiment, the calibrated extended image gather in the smearingfunction may be obtained in the slowness or angle domain with a Radontransform or source-direction gathers.

In an embodiment, the aggregate blurring function, a rate of change in ashape of the aggregate blurring function, the one or more plane-wavereflection coefficients, or a combination thereof may be displayed.

In an embodiment, the shot record may be generated by a shot fired by auser on land or in a marine environment.

Embodiments of the disclosure may provide a non-transitorycomputer-readable medium storing instructions that, when executed by atleast one processor of a computing system, cause the computing system toperform operations is disclosed. The operations may include determininga source wave field at a reference depth proximate to a reflector for ashot record. A receiver-side blurring functions may be determined at thereference depth. An aggregate blurring function may be constructed basedat least partially on the source wave field and the receiver-sideblurring function. A plane-wave reflection coefficient may be determinedbased at least partially on the aggregate blurring function.

Embodiments of the disclosure may provide a computing system including aprocessor and a memory system that includes a non-transitory,computer-readable medium including instructions that, when executed byat least one of the processor, cause the computing system to performoperations. The operations may include determining a source wave fieldat a reference depth proximate to a reflector for a shot record. Areceiver-side blurring function may be determined at the referencedepth. An aggregate blurring function may be constructed based at leastpartially on the source wave field and the receiver-side blurringfunction. A plane-wave reflection coefficient may be determined based atleast partially on the aggregate blurring function.

Embodiments of the disclosure may provide a computing system includingmeans for determining a source wave field at a reference depth proximateto a reflector for a shot record. The system also includes means fordetermining a receiver-side blurring function at the reference depth.The system also includes means for constructing an aggregate blurringfunction based at least partially on the source wave field and thereceiver-side blurring function. The system also includes means fordetermining a plane-wave reflection coefficient based at least partiallyon the aggregate blurring function.

Embodiments of the disclosure may provide a computing system including aprocessor and a memory system that includes a non-transitory,computer-readable medium including instructions that, when executed byat least one of the processor, cause the computing system to determine asource wave field at a reference depth proximate to a reflector for ashot record. The instructions also cause the system to determine areceiver-side blurring function at the reference depth. The instructionsalso cause the system to construct an aggregate blurring function basedat least partially on the source wave field and the receiver-sideblurring function. The instructions also cause the system to determine aplane-wave reflection coefficient based at least partially on theaggregate blurring function.

This summary is provided to introduce a selection of concepts that arefurther described below in the detailed description. This summary is notintended to identify key or essential features of the claimed subjectmatter, nor is it intended to be used as an aid in limiting the scope ofthe claimed subject matter.

BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings, which are incorporated in and constitute apart of this specification, illustrate embodiments of the presentteachings and together with the description, serve to explain theprinciples of the present teachings. In the figures:

FIGS. 1A, 1B, 1C, 1D, 2, 3A, and 3B illustrate simplified, schematicviews of an oilfield and its operation, according to an embodiment.

FIG. 4 illustrates a flowchart of a method for determining a reflectioncoefficient, according to an embodiment.

FIG. 5 illustrates another flowchart of a method for determining areflection coefficient, according to an embodiment.

FIG. 6 illustrates a two-dimensional view of an overburden and a targetinterface, according to an embodiment.

FIG. 7 illustrates a schematic view of a receiver-side blurringfunction, according to an embodiment.

FIG. 8 illustrates a two-dimensional view of a subterranean modelincluding first and second subsalt image points, according to anembodiment.

FIGS. 9 and 10 illustrate three-dimensional views of receiver-sideblurring functions centered on the first and second image points shownin FIG. 8, respectively, according to an embodiment.

FIGS. 11 and 12 illustrate three-dimensional views of Radon transformsor slant stacks of the receiver-side blurring functions shown in FIGS. 9and 10, respectively, according to an embodiment.

FIGS. 13 and 14 illustrate three-dimensional views of an aggregateblurring function corresponding to a circular salt body in an overburdenand a horizontal subsalt reflector, respectively, according to anembodiment.

FIG. 15 illustrates an extended image gather for a horizontal interfaceplotted in the (x,t) plane, according to an embodiment.

FIG. 16 illustrates a reflection operator for the interface shown inFIG. 15 plotted in the same manner, according to an embodiment.

FIG. 17 illustrates a synthetic extended image gather obtained byapplying the aggregate blurring function to the reflection operatorshown in FIG. 16, according to an embodiment.

FIGS. 18A-C illustrate a flowchart of a method for determining areflection coefficient, according to an embodiment.

FIG. 19 illustrates a computing system for implementing one or more ofthe methods disclosed herein, according to an embodiment.

DESCRIPTION OF EMBODIMENTS

Reference will now be made in detail to embodiments, examples of whichare illustrated in the accompanying drawings and figures. In thefollowing detailed description, numerous specific details are set forthin order to provide a thorough understanding of the invention. However,it will be apparent to one of ordinary skill in the art that theinvention may be practiced without these specific details. In otherinstances, well-known methods, procedures, components, circuits andnetworks have not been described in detail so as not to unnecessarilyobscure aspects of the embodiments.

It will also be understood that, although the terms first, second, etc.may be used herein to describe various elements, these elements shouldnot be limited by these terms. These terms are only used to distinguishone element from another. For example, a first object or step could betermed a second object or step, and, similarly, a second object or stepcould be termed a first object or step, without departing from the scopeof the invention. The first object or step, and the second object orstep, are both, objects or steps, respectively, but they are not to beconsidered the same object or step.

The terminology used in the description of the invention herein is forthe purpose of describing particular embodiments only and is notintended to be limiting of the invention. As used in the description ofthe invention and the appended claims, the singular forms “a,” “an” and“the” are intended to include the plural forms as well, unless thecontext clearly indicates otherwise. It will also be understood that theterm “and/or” as used herein refers to and encompasses any and allpossible combinations of one or more of the associated listed items. Itwill be further understood that the terms “includes,” “including,”“comprises” and/or “comprising,” when used in this specification,specify the presence of stated features, integers, steps, operations,elements, and/or components, but do not preclude the presence oraddition of one or more other features, integers, steps, operations,elements, components, and/or groups thereof. Further, as used herein,the term “if” may be construed to mean “when” or “upon” or “in responseto determining” or “in response to detecting,” depending on the context.

Attention is now directed to processing procedures, methods, techniquesand workflows that are in accordance with some embodiments. Someoperations in the processing procedures, methods, techniques andworkflows disclosed herein may be combined and/or the order of someoperations may be changed.

FIGS. 1A-1D illustrate simplified, schematic views of oilfield 100having subterranean formation 102 containing reservoir 104 therein inaccordance with implementations of various technologies and techniquesdescribed herein. FIG. 1A illustrates a survey operation being performedby a survey tool, such as seismic truck 106.1, to measure properties ofthe subterranean formation. The survey operation is a seismic surveyoperation for producing sound vibrations. In FIG. 1A, one such soundvibration, e.g., sound vibration 112 generated by source 110, reflectsoff horizons 114 in earth formation 116. A set of sound vibrations isreceived by sensors, such as geophone-receivers 118, situated on theearth's surface. The data received 120 is provided as input data to acomputer 122.1 of a seismic truck 106.1, and responsive to the inputdata, computer 122.1 generates seismic data output 124. This seismicdata output may be stored, transmitted or further processed as desired,for example, by data reduction.

FIG. 1B illustrates a drilling operation being performed by drillingtools 106.2 suspended by rig 128 and advanced into subterraneanformations 102 to form wellbore 136. Mud pit 130 is used to drawdrilling mud into the drilling tools via flow line 132 for circulatingdrilling mud down through the drilling tools, then up wellbore 136 andback to the surface. The drilling mud is typically filtered and returnedto the mud pit. A circulating system may be used for storing,controlling, or filtering the flowing drilling mud. The drilling toolsare advanced into subterranean formations 102 to reach reservoir 104.Each well may target one or more reservoirs. The drilling tools areadapted for measuring downhole properties using logging while drillingtools. The logging while drilling tools may also be adapted for takingcore sample 133 as shown.

Computer facilities may be positioned at various locations about theoilfield 100 (e.g., the surface unit 134) and/or at remote locations.Surface unit 134 may be used to communicate with the drilling toolsand/or offsite operations, as well as with other surface or downholesensors. Surface unit 134 is capable of communicating with the drillingtools to send commands to the drilling tools, and to receive datatherefrom. Surface unit 134 may also collect data generated during thedrilling operation and produce data output 135, which may then be storedor transmitted.

Sensors (S), such as gauges, may be positioned about oilfield 100 tocollect data relating to various oilfield operations as describedpreviously. As shown, sensor (S) is positioned in one or more locationsin the drilling tools and/or at rig 128 to measure drilling parameters,such as weight on bit, torque on bit, pressures, temperatures, flowrates, compositions, rotary speed, and/or other parameters of the fieldoperation. Sensors (S) may also be positioned in one or more locationsin the circulating system.

Drilling tools 106.2 may include a bottom hole assembly (BHA) (notshown), generally referenced, near the drill bit (e.g., within severaldrill collar lengths from the drill bit). The bottom hole assemblyincludes capabilities for measuring, processing, and storinginformation, as well as communicating with surface unit 134. The bottomhole assembly further includes drill collars for performing variousother measurement functions.

The bottom hole assembly may include a communication subassembly thatcommunicates with surface unit 134. The communication subassembly isadapted to send signals to and receive signals from the surface using acommunications channel such as mud pulse telemetry, electro-magnetictelemetry, or wired drill pipe communications. The communicationsubassembly may include, for example, a transmitter that generates asignal, such as an acoustic or electromagnetic signal, which isrepresentative of the measured drilling parameters. It will beappreciated by one of skill in the art that a variety of telemetrysystems may be employed, such as wired drill pipe, electromagnetic orother known telemetry systems.

Typically, the wellbore is drilled according to a drilling plan that isestablished prior to drilling. The drilling plan typically sets forthequipment, pressures, trajectories and/or other parameters that definethe drilling process for the wellsite. The drilling operation may thenbe performed according to the drilling plan. However, as information isgathered, the drilling operation may deviate from the drilling plan.Additionally, as drilling or other operations are performed, thesubsurface conditions may change. The earth model may also adjust as newinformation is collected

The data gathered by sensors (S) may be collected by surface unit 134and/or other data collection sources for analysis or other processing.The data collected by sensors (S) may be used alone or in combinationwith other data. The data may be collected in one or more databasesand/or transmitted on or offsite. The data may be historical data, realtime data, or combinations thereof. The real time data may be used inreal time, or stored for later use. The data may also be combined withhistorical data or other inputs for further analysis. The data may bestored in separate databases, or combined into a single database.

Surface unit 134 may include transceiver 137 to allow communicationsbetween surface unit 134 and various portions of the oilfield 100 orother locations. Surface unit 134 may also be provided with orfunctionally connected to one or more controllers (not shown) foractuating mechanisms at oilfield 100. Surface unit 134 may then sendcommand signals to oilfield 100 in response to data received. Surfaceunit 134 may receive commands via transceiver 137 or may itself executecommands to the controller. A processor may be provided to analyze thedata (locally or remotely), make the decisions and/or actuate thecontroller. In this manner, oilfield 100 may be selectively adjustedbased on the data collected. This technique may be used to optimize (orimprove) portions of the field operation, such as controlling drilling,weight on bit, pump rates, or other parameters. These adjustments may bemade automatically based on computer protocol, and/or manually by anoperator. In some cases, well plans may be adjusted to select optimum(or improved) operating conditions, or to avoid problems.

FIG. 1C illustrates a wireline operation being performed by wirelinetool 106.3 suspended by rig 128 and into wellbore 136 of FIG. 1B.Wireline tool 106.3 is adapted for deployment into wellbore 136 forgenerating well logs, performing downhole tests and/or collectingsamples. Wireline tool 106.3 may be used to provide another method andapparatus for performing a seismic survey operation. Wireline tool 106.3may, for example, have an explosive, radioactive, electrical, oracoustic energy source 144 that sends and/or receives electrical signalsto surrounding subterranean formations 102 and fluids therein.

Wireline tool 106.3 may be operatively connected to, for example,geophones 118 and a computer 122.1 of a seismic truck 106.1 of FIG. 1A.Wireline tool 106.3 may also provide data to surface unit 134. Surfaceunit 134 may collect data generated during the wireline operation andmay produce data output 135 that may be stored or transmitted. Wirelinetool 106.3 may be positioned at various depths in the wellbore 136 toprovide a survey or other information relating to the subterraneanformation 102.

Sensors (S), such as gauges, may be positioned about oilfield 100 tocollect data relating to various field operations as describedpreviously. As shown, sensor S is positioned in wireline tool 106.3 tomeasure downhole parameters which relate to, for example porosity,permeability, fluid composition and/or other parameters of the fieldoperation.

FIG. 1D illustrates a production operation being performed by productiontool 106.4 deployed from a production unit or Christmas tree 129 andinto completed wellbore 136 for drawing fluid from the downholereservoirs into surface facilities 142. The fluid flows from reservoir104 through perforations in the casing (not shown) and into productiontool 106.4 in wellbore 136 and to surface facilities 142 via gatheringnetwork 146.

Sensors (S), such as gauges, may be positioned about oilfield 100 tocollect data relating to various field operations as describedpreviously. As shown, the sensor (S) may be positioned in productiontool 106.4 or associated equipment, such as Christmas tree 129,gathering network 146, surface facility 142, and/or the productionfacility, to measure fluid parameters, such as fluid composition, flowrates, pressures, temperatures, and/or other parameters of theproduction operation.

Production may also include injection wells for added recovery. One ormore gathering facilities may be operatively connected to one or more ofthe wellsites for selectively collecting downhole fluids from thewellsite(s).

While FIGS. 1B-1D illustrate tools used to measure properties of anoilfield, it will be appreciated that the tools may be used inconnection with non-oilfield operations, such as gas fields, mines,aquifers, storage or other subterranean facilities. Also, while certaindata acquisition tools are depicted, it will be appreciated that variousmeasurement tools capable of sensing parameters, such as seismic two-waytravel time, density, resistivity, production rate, etc., of thesubterranean formation and/or its geological formations may be used.Various sensors (S) may be located at various positions along thewellbore and/or the monitoring tools to collect and/or monitor thedesired data. Other sources of data may also be provided from offsitelocations.

The field configurations of FIGS. 1A-1D are intended to provide a briefdescription of an example of a field usable with oilfield applicationframeworks. Part of, or the entirety, of oilfield 100 may be on land,water and/or sea. Also, while a single field measured at a singlelocation is depicted, oilfield applications may be utilized with anycombination of one or more oilfields, one or more processing facilitiesand one or more wellsites.

FIG. 2 illustrates a schematic view, partially in cross section ofoilfield 200 having data acquisition tools 202.1, 202.2, 202.3 and 202.4positioned at various locations along oilfield 200 for collecting dataof subterranean formation 204 in accordance with implementations ofvarious technologies and techniques described herein. Data acquisitiontools 202.1-202.4 may be the same as data acquisition tools 106.1-106.4of FIGS. 1A-1D, respectively, or others not depicted. As shown, dataacquisition tools 202.1-202.4 generate data plots or measurements208.1-208.4, respectively. These data plots are depicted along oilfield200 to demonstrate the data generated by the various operations.

Data plots 208.1-208.3 are examples of static data plots that may begenerated by data acquisition tools 202.1-202.3, respectively; however,it should be understood that data plots 208.1-208.3 may also be dataplots that are updated in real time. These measurements may be analyzedto better define the properties of the formation(s) and/or determine theaccuracy of the measurements and/or for checking for errors. The plotsof each of the respective measurements may be aligned and scaled forcomparison and verification of the properties.

Static data plot 208.1 is a seismic two-way response over a period oftime. Static plot 208.2 is core sample data measured from a core sampleof the formation 204. The core sample may be used to provide data, suchas a graph of the density, porosity, permeability, or some otherphysical property of the core sample over the length of the core. Testsfor density and viscosity may be performed on the fluids in the core atvarying pressures and temperatures. Static data plot 208.3 is a loggingtrace that typically provides a resistivity or other measurement of theformation at various depths.

A production decline curve or graph 208.4 is a dynamic data plot of thefluid flow rate over time. The production decline curve typicallyprovides the production rate as a function of time. As the fluid flowsthrough the wellbore, measurements are taken of fluid properties, suchas flow rates, pressures, composition, etc.

Other data may also be collected, such as historical data, user inputs,economic information, and/or other measurement data and other parametersof interest. As described below, the static and dynamic measurements maybe analyzed and used to generate models of the subterranean formation todetermine characteristics thereof. Similar measurements may also be usedto measure changes in formation aspects over time.

The subterranean structure 204 has a plurality of geological formations206.1-206.4. As shown, this structure has several formations or layers,including a shale layer 206.1, a carbonate layer 206.2, a shale layer206.3 and a sand layer 206.4. A fault 207 extends through the shalelayer 206.1 and the carbonate layer 206.2. The static data acquisitiontools are adapted to take measurements and detect characteristics of theformations.

While a specific subterranean formation with specific geologicalstructures is depicted, it will be appreciated that oilfield 200 maycontain a variety of geological structures and/or formations, sometimeshaving extreme complexity. In some locations, typically below the waterline, fluid may occupy pore spaces of the formations. Each of themeasurement devices may be used to measure properties of the formationsand/or its geological features. While each acquisition tool is shown asbeing in specific locations in oilfield 200, it will be appreciated thatone or more types of measurement may be taken at one or more locationsacross one or more fields or other locations for comparison and/oranalysis.

The data collected from various sources, such as the data acquisitiontools of FIG. 2, may then be processed and/or evaluated. Typically,seismic data displayed in static data plot 208.1 from data acquisitiontool 202.1 is used by a geophysicist to determine characteristics of thesubterranean formations and features. The core data shown in static plot208.2 and/or log data from well log 208.3 are typically used by ageologist to determine various characteristics of the subterraneanformation. The production data from graph 208.4 is typically used by thereservoir engineer to determine fluid flow reservoir characteristics.The data analyzed by the geologist, geophysicist and the reservoirengineer may be analyzed using modeling techniques.

FIG. 3A illustrates an oilfield 300 for performing production operationsin accordance with implementations of various technologies andtechniques described herein. As shown, the oilfield has a plurality ofwellsites 302 operatively connected to central processing facility 354.The oilfield configuration of FIG. 3A is not intended to limit the scopeof the oilfield application system. Part, or all, of the oilfield may beon land and/or sea. Also, while a single oilfield with a singleprocessing facility and a plurality of wellsites is depicted, anycombination of one or more oilfields, one or more processing facilitiesand one or more wellsites may be present.

Each wellsite 302 has equipment that forms wellbore 336 into the earth.The wellbores extend through subterranean formations 306 includingreservoirs 304. These reservoirs 304 contain fluids, such ashydrocarbons. The wellsites draw fluid from the reservoirs and pass themto the processing facilities via surface networks 344. The surfacenetworks 344 have tubing and control mechanisms for controlling the flowof fluids from the wellsite to processing facility 354.

Attention is now directed to FIG. 3B, which illustrates a side view of amarine-based survey 360 of a subterranean subsurface 362 in accordancewith one or more implementations of various techniques described herein.Subsurface 362 includes seafloor surface 364. Seismic sources 366 mayinclude marine sources such as vibroseis or airguns, which may propagateseismic waves 368 (e.g., energy signals) into the Earth over an extendedperiod of time or at a nearly instantaneous energy provided by impulsivesources. The seismic waves may be propagated by marine sources as afrequency sweep signal. For example, marine sources of the vibroseistype may initially emit a seismic wave at a low frequency (e.g., 5 Hz)and increase the seismic wave to a high frequency (e.g., 80-90 Hz) overtime.

The component(s) of the seismic waves 368 may be reflected and convertedby seafloor surface 364 (i.e., reflector), and seismic wave reflections370 may be received by a plurality of seismic receivers 372. Seismicreceivers 372 may be disposed on a plurality of streamers (i.e.,streamer array 374). The seismic receivers 372 may generate electricalsignals representative of the received seismic wave reflections 370. Theelectrical signals may be embedded with information regarding thesubsurface 362 and captured as a record of seismic data.

In one implementation, each streamer may include streamer steeringdevices such as a bird, a deflector, a tail buoy and the like, which arenot illustrated in this application. The streamer steering devices maybe used to control the position of the streamers in accordance with thetechniques described herein.

In one implementation, seismic wave reflections 370 may travel upwardand reach the water/air interface at the water surface 376, a portion ofreflections 370 may then reflect downward again (i.e., sea-surface ghostwaves 378) and be received by the plurality of seismic receivers 372.The sea-surface ghost waves 378 may be referred to as surface multiples.The point on the water surface 376 at which the wave is reflecteddownward is generally referred to as the downward reflection point.

The electrical signals may be transmitted to a vessel 380 viatransmission cables, wireless communication or the like. The vessel 380may then transmit the electrical signals to a data processing center. Inanother embodiment, the vessel 380 may include an onboard computercapable of processing the electrical signals (i.e., seismic data). Thoseskilled in the art having the benefit of this disclosure will appreciatethat this illustration is highly idealized. For instance, surveys may beof formations deep beneath the surface. The formations may typicallyinclude multiple reflectors, some of which may include dipping events,and may generate multiple reflections (including wave conversion) forreceipt by the seismic receivers 372. In one implementation, the seismicdata may be processed to generate a seismic image of the subsurface 362.

Typically, marine seismic acquisition systems tow each streamer instreamer array 374 at the same depth (e.g., 5-10 m). However, marinebased survey 360 may tow each streamer in streamer array 374 atdifferent depths such that seismic data may be acquired and processed ina manner that avoids the effects of destructive interference due tosea-surface ghost waves. For instance, marine-based survey 360 of FIG.3B illustrates eight streamers towed by vessel 380 at eight differentdepths. The depth of each streamer may be controlled and maintainedusing the birds disposed on each streamer.

FIG. 4 illustrates a flowchart of a method 400 for determining areflection coefficient, according to an embodiment. The method 400 maybegin by selecting an image point x′ and a reference depth z_(ref) in asubterranean formation, as at 402. The reference depth may be just abovea nearby reflector and/or may be at the actual reflector depth z_(R) ifit is a horizontal reflector. One or more source wave fields may bedetermined at the reference depth z_(ref) for one or more shot records,as at 404. In addition, one or more receiver-side blurring functions maybe determined at the reference depth z_(ref), as at 406. An aggregateblurring function may be constructed based on (e.g., by combining) theone or more source wave fields and the one or more receiver-sideblurring functions, as at 408. Once the last shot record has occurred,one or more reflection coefficients may be determined based on theaggregate blurring function, as at 410. During the process, theaggregate blurring function may be decomposed into plane-wavecomponents, and it may be defined and used over a range of depths belowthe reference depth z_(ref). Furthermore, one may use this function ineither the time or frequency domain.

FIG. 5 illustrates a flowchart of a method 500 for determining aplane-wave reflection coefficient, according to an embodiment that usesthe frequency-domain option for exploiting the aggregate blurringfunction. The method 500 may include selecting an image point x′ and areference depth z_(ref) in a subterranean formation, as at 502. One ormore source wave fields may be determined proximate to the image pointx′ at the reference depth z_(ref) for one or more shot records, as at504. The one or more shot records may be generated by one or more shotsfired by a user. The shots may be fired on land or at sea. The sourcewave field may be a downgoing source wave field in the space-time domain(x,t). The source wave field may be transformed into the frequencydomain, as at 506.

One or more receiver-side blurring functions may be determined proximateto the image point x′ at the reference depth z_(ref), as at 508. Thereceiver-side blurring functions may be in the space-time domain (x,t).The one or more receiver-side blurring functions may be transformed intothe frequency domain, as at 510.

An aggregate blurring function may be constructed based at least in parton (e.g., by combining) the one or more source wave fields and the oneor more receiver-side blurring functions, as at 512. The aggregateblurring function may be in the frequency domain if its components arein the frequency domain.

Migration data may be received, as at 514. In at least one embodiment,the migration data may be or include depth migration data obtained byreverse-time migration (“RTM”) or any method that permits the formationof the subsurface extended image gathers referred to below. An extendedimage gather (“EIG”) may be constructed at the reference depth z_(ref)based at least in part on the migration data, as at 516. The extendedimage gather may be in the space-time domain (x,t). The extended imagegather may be transformed into the frequency domain, as at 518.

The extended image gather may be calibrated for a band limit of the oneor more receiver-side blurring functions, as at 520. The extended imagegather may also be calibrated for a well-tie for source-side andreceiver-side modelling uncertainty.

Once the last shot record has occurred, one or more reflection operatorsmay be determined based at least in part on the aggregate blurringfunction and the calibrated extended image gather, as at 522. The one ormore reflection operators may be in the space-frequency domain in thisembodiment, e.g., if the extended image gather and aggregate blurringfunctions are in this domain.

The one or more reflection operators may be transformed into theslowness and/or angle domain to give plane-wave reflection coefficients,as at 524. The plane-wave reflection coefficients may be used as inputsto an AVOA analysis for material properties contrasts at the reflectinginterface. These results may then be shown on a display (e.g., amonitor). For example, the aggregate blurring function may be displayedat various places to ascertain how wide and/or stable it is. A width ofthe blurring function in space may be used to determine the minimumlength scale that may be resolved in space. Two small scatters separatedby less than this distance may not be distinguished. A rate of change ina shape of the aggregate blurring function may determine how well theslowness/angle in the EIG may be resolved, meaning the plane-wavereflection coefficients “in” the EIG. In another embodiment, the EIGitself may be displayed to assess how prominent the reflector dip maybe. The one or more plane-wave reflection coefficients may also bedisplayed to ascertain rock properties and the presence of hydrocarbonsproximate to the reflector.

This method 500 may account for the laws of reflection and compensatefor smearing of the plane-wave reflection coefficients in the slownessand/or angle domain as well as spatial-domain illumination and blurring.

Depth Imaging

Depth imaging may use a velocity model for the survey area. This modelmay be smoothly varying, but major refractive bodies such as salt may beexplicitly contained if reflecting interfaces below these bodies.

The depth-migration technique may provide gridded estimates of theincident and reflected waves near the target interface and of thereceiver-side blurring function. One such technique is calledreverse-time migration (“RTM”), and it uses wave fields sampled on agrid in the migration velocity model. The gridded wave fields near thereflector may be obtained, for example, by a finite-difference or aray-based extrapolation technique.

The incident-wave estimate at the interface may be obtained byextrapolation from the known source location in the migration velocitymodel. The reflected-wave estimate at the interface may be obtained byreverse-time extrapolation of the field recordings for the same shot,which may be fed into the model at the surface-receiver locations.

Extended Image Gather (“EIG”)

FIG. 6 illustrates a two-dimensional view of an overburden 602 and ahorizontal target interface 604, according to an embodiment. Theoverburden 602 and a locally homogeneous background 606 are shown near(e.g., above) the target interface 604. The dashed contours represent aspace-domain blurring function.

Spatial coordinates (x,z) are used, plus time t. In the vicinity of theinterface, the estimated incident wave is denoted by d(x,z,t), and theestimated reflected wave by ũ(x,z,t). The migration tool described abovemay be used to obtain the reflected-wave estimate ũ(x,z,t). The localform of the analytical back propagator defined below may also beinvolved.

The EIG denoted by I(x, x′, z, t) and formed from these wavefields maybe defined by

$\begin{matrix}{{{I\left( {x,x^{\prime},z,t} \right)} = {\sum\limits_{s}{{\overset{\sim}{u}\left( {x,z,{t - t^{\prime}}} \right)}{d^{\dagger}\left( {x^{\prime},z,t^{\prime}} \right)}}}},} & (1)\end{matrix}$where superscript † denotes a time-reversed field, and summation is overshot index s. Integration with respect to time may be implicit on theright-hand side. The repeated variables may be integrated. Hence, theEIG may be a temporal correlation of the wavefields at spatiallyseparated points x and x′ at a given depth z.

The point (x′, z) in the arguments of the EIG I(x, x′, z, t) is calledthe image point, and it may be regarded as a constant for most of thefollowing. The EIG I(x, x′, z, t) may appear in the (x,t) plane at fixedz and x′. This (x,t) section may be called the subsurface“virtual-survey,” because it relates closely to the reflection operatorkernel, and the kernel may resemble the response due to a particularsource placed just above the interface.

In another embodiment, I(x, x′, z, t) may be in the (x,z) plane at afixed output time, (e.g., t=0), and this embodiment is also consideredbelow.

Mathematical Model of the EIG

Considering a flat reflector as a simple example let z_(R) be the depthof the interface (see FIG. 6). The medium immediately above theinterface may be locally homogeneous.

The interface reflection operator R(x,x′, t) may connect an incidentwave d(x′, z_(R), t) and the true upward-going wave u(x, z_(R), t) atthe reflector according tou(x,z _(R) ,t)=R(x,x′,t−t′)d(x′,z _(R) ,t′),  (2)where there is implicit integration over x′ and t′. Note that x′ inEquation (2) is a variable of integration and not the image-pointargument x′ in the EIG I(x, x′, z, t). A single t is displayed for time(i.e., temporal convolutions are implicit, with the focus on spatialvariables); however, the symbol ‘•’ may be used instead of t, in orderto increase further the spatial emphasis and to allow forfrequency-domain representations.

A simplified mathematical expression for the space-time domainreflection operator at the interface is

$\begin{matrix}{{{R\left( {x,x^{\prime},t} \right)} = {{- \frac{\overset{.}{\overset{\_}{\delta}}(t)}{2\;\pi}}*{\int{{R(p)}{\delta\left( {t - {p\left( {x - x^{\prime}} \right)}} \right)}{\mathbb{d}p}}}}},} & (3)\end{matrix}$

where R(p) is a plane-wave reflection coefficient at slowness p, anoverbar denotes a Hilbert transform, and a dot a time derivative. Theplane-wave reflection coefficient R(p) may be considered as part of anAVOA analysis. Further, when estimating R(p) from an EIG, blurring inspace and slowness/angle may be accounted for, as will be describedbelow.

The availability of an accurate gridded incident wave d (x′, z, t) maybe assumed around the reflector for each shot. This may useforward-in-time modelling.

On the receiver side, the corresponding reflected wave u(x, z_(R), t) atthe interface may not be available. An estimate ũ(x,z,t) may beavailable, e.g., as obtained by back-propagation or reverse-timeextrapolation of the field-data surface recordings.

The relation between these two types of reflected wave may berepresented by{tilde over (u)}(x,z,t)=W ^(−(†))(x,z,x _(r) ,z _(r) ,t)a(x _(r))W ⁻(x_(r) ,z _(r) ,x″,z _(R) ,t)u(x″,z _(R) ,t),  (4)where W⁻(x_(r), z_(r), x″, z_(R), t) is the upgoing-wave propagator fromthe reflector to the receiver level, z_(r), a(x_(r)) is a receiveraperture function of position x_(r) at that level, and W^(−(†))(x, z,x_(r), z_(r), t) is the back propagator from the receiver level to theEIG level. The latter may be a reversed-wavepath operator rather thansimply the time-reversed form of W⁻, hence the superscript symbol (†).As explained below, the forward and backward wave propagators may beobtained by full wave equation modelling, and they contain overburdentransmission losses as well as geometrical spreading. These may besurface to surface wave propagators, and they may be equivalent totrue-amplitude transmission operators for the region between the targetzone and the receivers. These basic propagators may contain internalmultiples. Even more complicated wavepaths may be dealt with thatencompass energy travelling deeper than the target for some portion ofthe wavepath. However, considering, as an example, the primary or directwavepath from the surface to depth and vice versa, an appropriatespace-time waveform windowing may be implicit.

From Equation (4), a receiver-side blurring function may be defined byB _(r)(x,z,x″,z _(R) ,t)=W ^(−(†))(x,z,x _(r) ,z _(r) ,t)a(x _(r))W ⁻(x_(r) ,z _(r) ,x″,z _(R) ,t).  (5)B_(r) is a function of two depths as well as two lateral variables, andtime. Collecting together the above information, the model for EIGanalysis may be written

$\begin{matrix}{{{I\left( {x,x^{\prime},z,t} \right)} = {\sum\limits_{s}{{B_{r}\left( {x,z,x^{''},z_{R},t} \right)}{R\left( {x^{''},x^{\prime\prime\prime},t} \right)}{d\left( {x^{\prime\prime\prime},z_{R},t} \right)}{d^{\dagger}\left( {x^{\prime},z,t} \right)}}}},} & (6)\end{matrix}$where on the left hand side, I(x, x′, z, t) is a gridded data objectgiven by Equation (1), and on the right hand side, R(x, x′, t) is amathematical model of the reflection operator containing plane-wavereflection coefficients according to Equation (3). There may beintegration over the dummy variables x″ and x′″ and/or timeconvolutions.

Changing the orders of integration and summation in Equation (6) mayyieldI(x,x′,z,t)=

(x,x′,x″,x′″,z,z _(R) ,t)R(x″,x′″,t),  (7)where

(x,x′,x″,x′″,z,z _(R) ,t)=Σ_(s) B _(r)(x,z,x″,z _(R) ,t)d(x′″,z _(R),t)d ^(†)(x′,z,t),   (8)may be an “aggregate blurring function,” which may contain at least someof the effects of the shots. The aggregate blurring function may alsoaccount for at least some of the effects of source-side illuminationvariations, as well as blurring due to back propagation offinite-aperture receiver data. z_(R) may be implicit in the reflectionoperator, and there may be no integration over the repeated depthvariables. There may be integration over the repeated lateral spatialvariables and time. Also, the receiver-side blurring function B_(r) maydepend at least partially on the shot index s for a moving array such asa marine streamer.

The aggregate blurring function

(Equation (8)) may be considered a function of three lateral spacecoordinates: x, x″ and x′″, as x′ in the EIG I(x, x′, z, t) may beconsidered a constant (e.g., the image point at which the reflectivityis to be found). The first coordinate x is the running coordinate of theEIG, and the other two are connected to the reflection operator.

The EIG on the left hand side of Equation (6) or (7) contains the sourcesignature twice, once in the estimated incident wave d^(†)(x′, z, t) andonce in the back-propagated field data ũ(x,z,t). On the right hand sideof these equations, the source signature again appears twice, in thesource-side combination d(x′″, z_(R), t)d^(†)(x′, z, t).

Thus, Equations (6) and (7) may be balanced with respect to thebandwidth of the source. At this point, the receiver-side blurringfunction B_(r) (x, z, x″, z_(R), t) may have infinite bandwidth,notwithstanding its use of a back propagator with a particular smoothingeffect on evanescent waves. B_(r) (x, z, x″, z_(R), t) may be thought ofas being singular like the generalized function that is the reflectionoperator in Equation (3). This operator may contain illuminationeffects, but as yet it may not have a frequency limitation.

Numerical Estimation of the Receiver-Side Blurring Function and theCalibration Field

FIG. 7 illustrates a schematic view of a receiver-side blurring function700, according to an embodiment. The receiver-side blurring functionB_(r)(x, z, x″, z_(R), t) of Equation (5) may be numerically estimatedusing the modelling engine for the given migration velocity model andacquisition geometry parameters. This may be done by injecting a smallsource 702 around or near the image point, so that x″ may now be takento characterize the location of this injected source (note that x″ isdenoted x′ in FIG. 7). This injected source 702 may be bandlimited, andit may be controlled by the user independently of the field data in theEIG.

Another consideration is that the wave propagators in the receiver-sideblurring function Equation (5) are surface-to-surface field-to-fieldoperators. Referring to FIG. 7, let z_(ref) be a reference depth orlevel just above the support of the finite source injection volume. Letthe resulting injected (e.g., upgoing) wavefield 704 at z_(ref) bedenoted b(x, x′, z_(ref),t). This plays the role of a calibration field.

The injected wavefield 704 may be modeled up to a receiver array 706 andthen reinjected as a downgoing wave 708. This downgoing wave 708measured at the same level z_(ref) may be an estimate of the smoothedblurring function for level z_(ref) defined by{hacek over (B)} _(r)(x,z _(ref) ,x′,z _(ref) ,t)=B _(r)(x,z _(ref),x″,z _(ref) ,t)b(x″,x′,z _(ref) ,t),   (9)where a temporal convolution is also understood. This smoothed blurringfunction may be gridded and may be available for various values of x′.Further, x′ is a local variable in this equation, and is not to beconfused with the image-point argument x′ in the EIG.

Equation (9) involves integration with respect to x″, the second spatialargument of B_(r). This argument of B_(r) may be composed with the firstargument of the reflection operator in the EIG Equation (6) or (7). Byexploiting a reciprocity property, it is possible to transfer thesmoothing integration to the first argument of B_(r). This facilitatesthe combination of the field-data-independent smoothing expressed byEquation (9) with the EIG formula Equation (6) or (7). This may bereferred to as “matching of the bandwidths.”

Reciprocity of the Receiver-Side Blurring Function

If the receiver-side blurring function had the following simpleoffset-only dependenceB _(r)(x,z _(ref) ,x′,z _(ref) ,t)=B _(r)(x−x′,z _(ref),0,z _(ref),t),   (10)then composition over x′ would be equivalent to composition over x by astraightforward change of integration variables. In such a case, theposition of b in Equation (9) may be moved to “the left of” B_(r), andcomposition of the latter with R in Equation (6) or (7) may proceed asbefore.

Some numerical examples, e.g., in the SIGSBEE model, indicate thatoffset-only dependence of the receiver-side blurring function may not berelied upon, although in others it may be. FIG. 8 illustrates atwo-dimensional view of a subterranean model 800 including first (e.g.,left) and second (e.g., right) subsalt image points 802, 804, accordingto an embodiment. FIGS. 9 and 10 illustrate three-dimensional views ofreceiver-side blurring functions 900, 1000 centered on the first andsecond image points 802, 804 shown in FIG. 8, respectively, according toan embodiment. As shown, the image points 802, 804 are separated by 7500feet, which is greater than the width of the well-focusedblurring-function spike 1002 in FIG. 10. The blurring function 1000 maybe stable with respect to changes in its image point 804 location thatare comparable to the width of its spike 1002, which is comparable tothe radius of the point 804 in FIG. 8. In contrast, for the image point802 on the left in FIG. 8, the spike 902 in the receiver-side blurringfunction 900 in FIG. 9 is not well-focused, and it changes its form whenits image point 802 is perturbed by a distance comparable to the size ofthe point 802 (i.e., it is unstable).

FIGS. 11 and 12 illustrate three-dimensional views of Radon transformsor slant stacks 1100, 1200 of the receiver-side blurring functions 900,1000 shown in FIGS. 9 and 10, respectively, according to an embodiment.A blurring function may be considered stable if it is approximatelyoffset-only dependent over a region comparable to the width of itscentral spike.

More generally, if the receiver-side blurring function had thereciprocity propertyB _(r)(x,z _(ref) ,x′,z _(ref) ,t)=B _(r)(x′,z _(ref) ,x,z _(ref) ,t),  (11)then the repositioning of b in Equation (9) may again be achieved.Further, there may be a reciprocity property, of the form:B _(r)(x,z _(ref) ,x′,z _(ref),•)=

^(†)(x,x″,•)B _(r) ^([†]†)(x′″,z _(ref) ,x″,z _(ref),•)

⁻¹(x′″,x′,•),  (12)where

is an operator related to the impedance operator for the locallyhomeogeneous medium around the target. Focusing on the spatial argumentsin Equation (12), it may be seen that on the right hand side, x isassociated via

^(†) with x″ in B_(r) ^([†]†) (i.e., the second spatial argument of thelatter). Similarly, on the right hand side, x′ is associated via

with x′″ in B_(r) ^([†]†) (i.e. the first spatial argument of thelatter). Thus, there is a reciprocity property. The complicatedsuperscripts on B_(r) ^([†]†) indicate time-reversal effects.Considering the blurring function for two different depths, rather thanhaving z_(ref) appear twice, then it may be apparent why the outertime-reversal superscript † appears. The other superscript symbol [†] isassociated with wavefield up/down splitting operators at the receivers,the details of which are omitted here.

The impedance operator and the related operator

may be defined in the locally homogeneous medium. Therefore, theseoperators may, in an embodiment, depend on the offset, but in others maydepend on the offset in addition to other factors.

Local Propagators and the Incident-Wave Composite Function

The fields at different z levels around the target may be related by thewave propagator for a locally homogeneous acoustic medium. Ignoring thelocal evanescent waves, these propagators may be represented bysimplified plane-wave integrals such as:

$\begin{matrix}{{W^{\pm}\left( {x,z,x^{\prime},z^{\prime},t} \right)} = {{- \frac{\overset{.}{\overset{\_}{\delta}}(t)}{2\;\pi}}*{\int{{\delta\left( {t - \left( {{\pm {q\left( {z - z^{\prime}} \right)}} + {p\left( {x - x^{\prime}} \right)}} \right)} \right)}{\mathbb{d}p}}}}} & (13)\end{matrix}$with real vertical slowness q≧0 and z, z′ chosen as appropriate fordown-/up-going waves. There is a similar definition for the backpropagators W^(±(†)).

An incident-wave composite field at the imaging depth z may be definedto beD(x,x′,z,•)=d(x,z,•)d ^(†)(x′,z,•),  (14)and then form the dual-depth combinationd(x,z _(R),•)d ^(†)(x′,z,•)=W ⁺(x,z _(R) ,x″,z _(ref),•))W ^(+†)(x″,z_(ref) ,x′″,z,•)D(x′″,x′,z,•).  (15)

Similarly, the receiver-side blurring function purely at z_(ref) and thedual-depth version for (z, z_(R)) may be related according to:B _(r)(x,z,x′,z _(R),•)=W ^(+†)(x,z,x″,z _(ref),)B _(r)(x″,z _(ref),x′″,z _(ref),•)×W ⁻(x′″,z _(ref) ,x′,z _(R),•).  (16)

In Equations (15) and (16), for scalar fields, the operators may berearranged without the matrix transpositions. The order of integrationmay simply be changed. The local wave propagators may be functions ofoffset-only, like the local impedance operator and normal operator.

Bandwidth Matching and the Calibrated EIG

The calibrated EIG may be defined as the composition of the calibrationfunction and the EIG according to:{hacek over (I)}(x,x′,z,•)=b ^(†)(x″,x,•)I(x″,x′,z,•),  (17)where the breve on {hacek over (I)} indicates the extra smoothing due tothe inclusion of the bandwidth of b. This may be an example of bandwidthmatching, as it establishes a commonality between the field data and themodelled receiver-side blurring function. The definition in Equation(17) involves integration over the first argument of b. The EIG functionI is given by Equation (6) and B_(r) by Equation (16). Noting the factthat b is also a function of offset-only in the locally homogeneousmedium, it is possible to interchange the orders of b and W^(+†).

Next, as B_(r) has the reciprocity property in Equation (12) and as thenormal operator

is also a function of offset-only, it is possible to interchange theorders of b and

^(†). This brings b into composition with B_(r) ^([†]†) (x′″, z_(ref),x″, z_(ref),•) shown in Equation (12) and, as noted there, it is theintegration with respect to the second argument x″ which arises in thatcomposition.

This composition may correspond to the definition of the smooth function{hacek over (B)} _(r) ^([†]†)(x′″,z _(ref) ,x′,z _(ref),•)=B _(r)^([†]†)(x′″,z _(ref) ,x″,z _(ref),•)b ^(†)(x″,x′,z _(ref),•),  (18)which is the time-reversed form of the band-limited receiver-sideblurring function at level z_(ref) that is available from the numericalmodelling. Thus, reciprocity of the blurring function may facilitateembedding the available smoothed function in Equation (18) into the EIGin Equation (6).

These calculations may make use of the operator plane-wave expansions,such as Equations (3) and (13), and a similar one for the impedanceoperator. These may include forward and backward Radon transforms.

A result of the calibrated EIG may be written in the explicit form:

$\begin{matrix}{{{\overset{\Cup}{I}\left( {x,x^{\prime},z,t} \right)} = {\frac{{\overset{.}{\overset{\_}{\delta}}(t)}*{\overset{.}{\overset{\_}{\delta}}(t)}}{\left( {2\;\pi} \right)^{2}}*{\int{\int{{R\left( p^{\prime} \right)}{i(p)}{i^{- 1}\left( p^{\prime} \right)} \times {\sum\limits_{s}{{{\overset{\sim}{\overset{\sim}{\overset{\Cup}{B}}}}^{{\lbrack\dagger\rbrack}\dagger}\left( {{- p^{\prime}},z_{ref},{- p},z_{ref},{t + {q^{\prime}\left( {z - z_{R}} \right)} - {q^{\prime}\left( {z_{R} - z_{ref}} \right)} + {q\left( {z - z_{ref}} \right)} + {px}}} \right)}*{\overset{\sim}{D}\left( {p^{\prime},z,x^{\prime},t} \right)}{\mathbb{d}p}{\mathbb{d}p^{\prime}}}}}}}}},} & (19)\end{matrix}$where i(p)=ρ/2q is half the plane-wave impedance. A tilde in Equation(19) represents a Radon transform, which for the receiver-side blurringfunction, may be applied to both of its lateral arguments whereas forthe incident-wave function it applies once. The fact that the sameslowness p′ appears in both functions is an expression of Snell's lawand may be enforced by the form of the reflection operator.

The smoothed EIG of Equation (19) may be balanced with respect to thesource wavelet of the field data, which appears twice on the far leftand the far right. It may also be balanced with respect to the bandwidthor wavelet of the numerically estimated receiver-side blurring function,which appears once on the left and once on the right (in {hacek over(B)}_(r) ^([†]†)). Thus, this formula provides a consistent basis forAVOA.

Slant Stacking and the Plane-Wave Reflection Coefficient SmearingFunction

As the p integral in Equation (19) is in the form of an inverse Radontransform, it may be seen that the slant stack of the EIG is in theform:

$\begin{matrix}{\mspace{79mu}{{{\overset{\Cup}{I}\left( {p,x^{\prime},z,t} \right)} = {\int{{S\left( {p,p^{\prime},x^{\prime},z,z_{R},t} \right)}{R\left( p^{\prime} \right)}{\mathbb{d}p^{\prime}}}}}\mspace{79mu}{where}}} & (20) \\{{S\left( {p,p^{\prime},x^{\prime},z,z_{R},t} \right)} = {{- {i(p)}}{i^{- 1}\left( p^{\prime} \right)}\frac{\overset{.}{\overset{\_}{\delta}}(t)}{2\;\pi}*{\sum\limits_{s}{{{\overset{\sim}{\overset{\sim}{\overset{\Cup}{B}}}}^{{\lbrack\dagger\rbrack}\dagger}\left( {{- p^{\prime}},z_{ref},{- p},z_{ref},{t + {q^{\prime}\left( {z - z_{R}} \right)} - {q^{\prime}\left( {z_{R} - z_{ref}} \right)} + {q\left( {z - z_{ref}} \right)}}} \right)}*{\overset{\sim}{D}\left( {p^{\prime},z,x^{\prime},t} \right)}}}}} & (21)\end{matrix}$is the plane-wave reflection coefficient smearing function in theslowness domain. This smearing function may be obtained from theaggregate blurring function and once found may be used to obtain theplane-wave reflection coefficient from the slant stack of the calibratedEIG via Equation (20). This is a type of deconvolution problem.

With a laterally homogeneous overburden and infinite shot and receiverapertures, the plane-wave reflection coefficient smearing function S inEquation (21) becomes proportional to the delta function δ(p−p′), andthe slowness integral in Equation (20) collapses to a single reflectioncoefficient at slowness p.

A laterally homogeneous overburden and infinite shot and receiverapertures lead to receiver-side and aggregate blurring functions withoffset-only dependence. These blurring functions are therefore preciselystable in the sense that their shapes do not change as the image pointchanges. The blurring over slowness/angle in Equation (20) arises fromthe double spatial dependence or observable instability of the aggregateblurring function, which leads to the displayed double slownessdependence on p and p′.

Corresponding Methods for the Offset-Depth (x,z) Domain at Fixed Timet=0

In order to treat (x,z) gathers at constant time, the z dependence maybe consolidated. This means that in the mathematical model of the EIGthe incident wave at a general depth is expressed in terms of its valueat z_(ref), which is then modified by explicit local wave propagators.The effect may be represented as:

$\begin{matrix}{{{\overset{\Cup}{I}\left( {x,x^{\prime},z,t} \right)} = {{- \frac{{\overset{.}{\overset{\_}{\delta}}(t)}*{\overset{.}{\overset{\_}{\delta}}(t)}*{\overset{.}{\overset{\_}{\delta}}(t)}}{\left( {2\;\pi} \right)^{3}}}*{\int{\int{\int{{R\left( p^{\prime} \right)}{i(p)}{i^{- 1}\left( p^{\prime} \right)} \times {\sum\limits_{s}{{{\overset{\sim}{\overset{\sim}{\overset{\Cup}{B}}}}^{{\lbrack\dagger\rbrack}\dagger}\left( {{- p^{\prime}},z_{ref},{- p},z_{ref},{t + {\left( {q + q^{''}} \right)\left( {z - z_{ref}} \right)} - {2{q^{\prime}\left( {z_{R} - z_{ref}} \right)}} + {px}}} \right)}*{{\overset{\sim}{\overset{\sim}{D}}}_{2}\left( {p^{\prime},p^{''},{t - {p^{''}x^{\prime}}}} \right)}{\mathbb{d}p}{\mathbb{d}p^{\prime}}{\mathbb{d}p^{''}}}}}}}}}},} & (22)\end{matrix}$where D₂ represents the incident wave combination previously denoted byD as it appears at fixed level z_(ref). Compared to Equation (19), thereis now an extra slowness integral, which represents the localpropagation in depth of the incident wave away from z_(ref) via theadditional vertical slowness that has been incorporated into thereceiver-side blurring function in order to consolidate the zdependence. Note that the left-hand side of Equation (22) is still thecalibrated data EIG obtained after depth migration, now evaluated atvarious depths for fixed time.

The slant stack of the (x,z) gather at fixed time may now be stated as

$\begin{matrix}{\mspace{85mu}{{{\overset{\Cup}{I}\left( {p,x^{\prime},z,t} \right)} = {\int{{S_{2}\left( {p,p^{\prime},x^{\prime},z,z_{R},t} \right)}{R\left( p^{\prime} \right)}{\mathbb{d}p^{\prime}}}}}\mspace{79mu}{where}}} & (23) \\{{S_{2}\left( {p,p^{\prime},x^{\prime},z,z_{R},t} \right)} = {\frac{{\overset{.}{\overset{\_}{\delta}}(t)}*{\overset{.}{\overset{\_}{\delta}}(t)}}{\left( {2\;\pi} \right)^{2}}*{\int{\frac{{i\left( {p\left( {\mu,p^{\prime}} \right)} \right)}{i^{- 1}\left( p^{\prime} \right)}}{\left( {q + q^{''}} \right){\partial_{p}\mu}}*{\sum\limits_{s}{{{\overset{\sim}{\overset{\sim}{\overset{\Cup}{B}}}}^{{\lbrack\dagger\rbrack}\dagger}\left( {{- p^{\prime}},z_{ref},{- {p\left( {\mu,p^{''}} \right)}},z_{ref},{t - {2\;{q^{\prime}\left( {z_{R} - z_{ref}} \right)}} - {\left( {q + q^{''}} \right)z_{ref}}}} \right)}*{{\overset{\sim}{\overset{\sim}{D}}}_{2}\left( {p^{\prime},p^{''},x^{\prime},{t - {p^{''}x^{\prime}}}} \right)}{\mathbb{d}p^{''}}}}}}}} & (24)\end{matrix}$is a second type of plane-wave reflection coefficient smearing function.This second type of smearing function is also obtained from theaggregate blurring function and once found may be used to obtain theplane-wave reflection coefficient from the slant stack of the calibratedEIG via Equation (23). This again is a type of deconvolution problem.

Implicit Compensation for Blurring in AVOA

As Equations (19) and (20) are linear in the plane-wave reflectioncoefficient, the effect of slowness/angle smearing may be incorporateddirectly into a method of AVOA analysis as an extension of thepoint-spread function approach. Such an inversion method works directlywith the blurred functions and may not use an explicit plane-wavereflection coefficient from which the slowness/angle-domain smearing hasfirst been removed or deconvolved. Thus, in some embodiments, notheoretical obstacles exist to actual AVOA inversion for the materialproperty contrasts from the calibrated EIGs using the blurring formulaspresented above. It remains an option to remove the blurring/smearingeffect from the calibrated EIG to obtain the explicit reflectionoperator or reflection coefficients, as described for a particularworkflow below.

Explicit Compensation for Blurring in AVOA

In some embodiments, the methods disclosed herein may obtain deblurredreflection coefficients as a function of incidence angle from an EIG, inparticular an EIG corresponding to a horizontal reflector. The methodsmay account for the blurring of reflection coefficients due to practicallimitations encountered during acquisition and imperfect illumination.

The reflection coefficients may be obtained via Equation (6). Thisequation expresses an (x,t) EIG, evaluated at or near the reflectordepth, in terms of the reflection operator for the interface. The twoobjects are connected by the aggregate blurring function

of Equation (8).

When R corresponds to a horizontal reflector embedded in a simplebackground, it may be regarded as a function of Δ=x″−x′″ (offset) inEquation (7), and the dimensionality of the problem is reduced. Using Δin Equation (6) or (7) considered at the reflector depth givesI(x,x′,z,t)=

(x,x′,Δ+x′″,x′″,z _(R) ,z _(R) ,t)R(Δ,t),  (22)where there is integration over x′″ as Equations (6) and (7) show. Thus,for a horizontal interface the aggregate blurring function

has in effect two lateral space dimensions (x, Δ), as x′ is consideredto be a constant (e.g., the image point at which the reflectivity is tobe found). Examples relating to Equation (22) follow in FIGS. 13-17.

FIG. 13 illustrates a three-dimensional view of an aggregate blurringfunction 1300 corresponding to a circular salt body in an overburden,and FIG. 14 illustrates a three-dimensional view of an aggregateblurring function 1400 corresponding to a circular salt body in ahorizontal subsalt reflector, according to an embodiment. The view inFIG. 13 shows the (x,t) plane, and the view in FIG. 14 shows the (Δ, t)plane.

FIG. 15 illustrates an extended image gather 1500 for a horizontalinterface plotted in the (x,t) plane, according to an embodiment. FIG.16 illustrates a reflection operator 1600 for the interface shown inFIG. 15 plotted in the same manner, according to an embodiment. FIG. 17illustrates a synthetic extended image gather 1700 obtained by applyingthe aggregate blurring function to the reflection operator 1600 shown inFIG. 16, according to an embodiment.

Exploiting the separability of the problem in the frequency domain mayreduce the dimensionality further. In the frequency domain, Equation(22) may be replaced by a matrix form in which the rows and columns ofthe matrix form of

are indexed by (x, Δ) samples. Solving these discrete equations for R(Δ,ω) and then discrete Fourier transforming with respect to Δ gives thefrequency-spatial wavenumber domain form R(k, ω). It may then bepossible to extract radial lines from the (k, ω,) plane, obtainingplane-wave reflection coefficients for slowness p. These plane-wavereflection coefficients may then be taken as input to AVOA analysis. Itremains an option to perform AVOA without explicitly removing theblurring function to obtain the deblurred reflection operator orreflection coefficients, as described in the previous section onimplicit compensation for blurring.

The reflector may be dipping with respect to the grid on which thewavefields are sampled and the EIG is defined. The aggregate blurringfunction may be unchanged, but the reflection operator in Equation (7)no longer simplifies from R(x″, x′″, t) to R(Δ, t). Thedipping-interface reflection operator still embodies Snell's law and,while an additional slant stack of the aggregate blurring functionbecomes involved, the smearing of plane-wave reflection coefficientsover slowness/angle may still be quantified using the methodology. Oncemore the method may be applied in the (x,t) at constant depth or the(x,z) plane at fixed t, and the option exists to perform operations inthe frequency domain.

The aggregate blurring function may be controlled by the overburdenstructure and acquisition parameters, whereas the reflection operatormay be local to the target reflector. The spatial variability of theaggregate blurring function may define the extent of blurring overslowness/angle of plane-wave reflection coefficients in the reflectionoperator. The terms stable and unstable may be used to characterize thenature of the aggregate blurring function, with the latter implying moreslowness/angle blurring.

The modelling may use a finite-difference method or a ray-based method.As long as the EIG is calibrated for the smoothing inherent in thereceiver-side blurring function, the manner in which the EIG istransformed into the slowness/angle domain may vary. Methods such asRadon transform, Fourier transform, or source-direction gathers may beused in principle.

The method is not limited to a two-dimensional (x,z) world. Theaggregate blurring function and the reflection operator may be definedsimilarly in a three-dimensional world, and their relationship remainsunchanged.

The decision where to fire the shots and/or place the receivers for dataacquisition may be aided by computing approximate aggregate blurringfunctions and slowness/angle domain smearing functions for anapproximate trial model of the subsurface constructed before the fieldwork. These functions may be used in the process of seismic surveydesign, where the aim is to find shot and receiver positions thatoptimize the subsurface illumination. This may mean optimizing extendedimage gathers for AVOA, but it may also be simply interpreted to meanoptimal illumination for the final stacked image. The method forquantifying slowness/angle domain smearing brings a new tool to surveydesign. For example, source-side wavefields and receiver-side blurringfunctions may be displayed to quality control the survey geometry andillumination.

Attention is now directed to FIGS. 18A-C, which are flow diagramsillustrating a method 1800 for determining a reflection coefficient,according to an embodiment. Some operations in the method 1800 may becombined and/or the order of some operations may be changed. Further,some operations in the method 1800 may be combined with aspects of theexample workflows of FIGS. 4 and 5, and/or the order of some operationsin methods 400 and/or 500 may be changed to account for incorporation ofaspects of the workflow illustrated by one or more of FIGS. 4 and 5.

The method 1800 may include determining one or more source wave fieldsat a reference depth proximate to a reflector for one or more shotrecords, as at 1802 (e.g., FIG. 5, 504; one or more source wave fieldsare determined). In an embodiment, the one or more shot records aregenerated by one or more shots fired by a user on land or in a marineenvironment, as at 1804 (e.g., FIG. 5, 504; the shot records aregenerated by one or more shots fired by a user).

The method 1800 may also include determining one or more receiver-sideblurring functions at the reference depth, as at 1806 (e.g., FIG. 5,508; one or more receiver-side blurring functions are determinedproximate to the image point at the reference depth). The method 1800may also include constructing an aggregate blurring function based atleast partially on the one or more source wave fields and the one ormore receiver-side blurring functions, as at 1808 (e.g., FIG. 5; 512; anaggregate blurring function is constructed based at least partially onthe one or more source wave fields and the one or more receiver-sideblurring functions).

The method 1800 may also include constructing a calibrated extendedimage gather proximate to the reference depth based at least partiallyon migration data and a calibration field including a band limit of theone or more receiver-side blurring functions, as at 1810 (e.g., FIG. 5,516; an extended image gather is constructed based at least partially onthe migration data). The reflector may be a dipping reflector withrespect to a coordinate frame of the one or more source wave fields andthe calibrated extended image gather, as at 1812 (e.g., FIG. 5, 516; thereflector may be a dipping reflector).

The method 1800 may also include transforming the calibrated extendedimage gather and the aggregate blurring function to a space-frequencydomain at a fixed depth, as at 1814 (e.g., FIG. 5, 518; the extendedimage gather may be transformed into the frequency domain). The method1800 may also include obtaining a reflection operator in thespace-frequency domain by matrix inversion, as at 1816 (e.g., FIG. 5,518; a reflection operator may be obtained in the space-frequencydomain).

The method 1800 may also include determining one or more plane-wavereflection coefficients based at least partially on the aggregateblurring function, as at 1818 (e.g., FIG. 5, 522; one or more reflectionoperators may be determined). The one or more plane-wave reflectioncoefficients may be at least partially based on an interaction of theaggregate blurring function with a reflection operator containing theplane-wave reflection coefficients, as at 1820 (e.g., FIG. 5, 522; theone or more plane-wave reflection coefficients may be at least partiallybased on an interaction of the aggregate blurring function with areflection operator containing the plane-wave reflection coefficients).The one or more plane-wave reflection coefficients may be determinedbased at least partially on the aggregate blurring function in aspace-time domain at a fixed depth proximate to the reference depth, asat 1822 (e.g., FIG. 5, 522; the one or more plane-wave reflectioncoefficients may be determined based at least partially on the aggregateblurring function in a space-time domain). The one or more plane-wavereflection coefficients may be determined in a domain of lateralposition and depth at a fixed time, as at 1824 (e.g., FIG. 5, 522; theone or more plane-wave reflection coefficients may be determined in adomain of lateral position and depth at a fixed time).

The method 1800 may also include transforming the reflection operatorfrom the space-frequency domain into the one or more plane-wavereflection coefficients by a Fourier transform, as at 1826 (e.g., FIG.5, 522; the reflection operator may be transformed from thespace-frequency domain into the one or more plane-wave reflectioncoefficients). The method 1800 may also include transforming theaggregate blurring function into a smearing function in a slowness orangle domain of the one or more plane-wave reflection coefficients, asat 1828 (e.g., FIG. 5, 524; the one or more reflection operators may betransformed into a slowness and/or angle domain). The smearing functionmay be determined in a space-time domain at a fixed depth proximate tothe reference depth, as at 1830 (e.g., FIG. 5, 524; the smearingfunction may be determined in a space-time domain at a fixed depthproximate to the reference depth). The smearing function may bedetermined in a domain of lateral position and depth at a fixed time, asat 1832 (e.g., FIG. 5, 524; the smearing function may be determined in adomain of lateral position and depth at a fixed time). The calibratedextended image gather in the smearing function may be obtained based atleast partially on a Radon transform or source-direction gathers, as at1834 (e.g., FIG. 5, 524; the calibrated extended image gather in thesmearing function may be obtained based at least partially on a Radontransform or source-direction gathers).

The method 1800 may also include transforming the calibrated extendedimage gather into the slowness or angle domain, as at 1836 (e.g., FIG.5, 524; the one or more reflection operators may be transformed into theslowness and/or angle domain). The method 1800 may also includedetermining the one or more plane-wave reflection coefficients based atleast partially on the smearing function in the slowness or angledomain, as at 1838 (e.g., FIG. 5, 524; the one or more plane-wavereflection coefficients may be determined based at least partially onthe smearing function in the slowness or angle domain). The method 1800may also include displaying the aggregate blurring function, a rate ofchange in a shape of the aggregate blurring function, the one or moreplane-wave reflection coefficients, or a combination thereof, as at 1840(e.g., FIG. 5, 524; displaying a portion of the method 1800).

In some embodiments, the methods 400, 500, and 1800 may be executed by acomputing system. FIG. 19 illustrates an example of such a computingsystem 1900, in accordance with some embodiments. The computing system1900 may include a computer or computer system 1901A, which may be anindividual computer system 1901A or an arrangement of distributedcomputer systems. The computer system 1901A includes one or moreanalysis modules 1902 that are configured to perform various tasksaccording to some embodiments, such as one or more methods disclosedherein (e.g., methods 400, 500, 1800, and/or combinations and/orvariations thereof). To perform these various tasks, the analysis module1902 executes independently, or in coordination with, one or moreprocessors 1904, which is (or are) connected to one or more storagemedia 1906. The processor(s) 1904 is (or are) also connected to anetwork interface 1907 to allow the computer system 1901A to communicateover a data network 1909 with one or more additional computer systemsand/or computing systems, such as 1901B, 1901C, and/or 1901D (note thatcomputer systems 1901B, 1901C and/or 1901D may or may not share the samearchitecture as computer system 1901A, and may be located in differentphysical locations, e.g., computer systems 1901A and 1901B may belocated in a processing facility, while in communication with one ormore computer systems such as 1901C and/or 1901D that are located in oneor more data centers, and/or located in varying countries on differentcontinents).

A processor can include a microprocessor, microcontroller, processormodule or subsystem, programmable integrated circuit, programmable gatearray, or another control or computing device.

The storage media 1906 can be implemented as one or morecomputer-readable or machine-readable storage media. Note that while inthe example embodiment of FIG. 19 storage media 1906 is depicted aswithin computer system 1901A, in some embodiments, storage media 1906may be distributed within and/or across multiple internal and/orexternal enclosures of computing system 1901A and/or additionalcomputing systems. Storage media 1906 may include one or more differentforms of memory including semiconductor memory devices such as dynamicor static random access memories (DRAMs or SRAMs), erasable andprogrammable read-only memories (EPROMs), electrically erasable andprogrammable read-only memories (EEPROMs) and flash memories, magneticdisks such as fixed, floppy and removable disks, other magnetic mediaincluding tape, optical media such as compact disks (CDs) or digitalvideo disks (DVDs), BLUERAY® disks, or other types of optical storage,or other types of storage devices. Note that the instructions discussedabove can be provided on one computer-readable or machine-readablestorage medium, or in another embodiment, can be provided on multiplecomputer-readable or machine-readable storage media distributed in alarge system having possibly plural nodes. Such computer-readable ormachine-readable storage medium or media is (are) considered to be partof an article (or article of manufacture). An article or article ofmanufacture can refer to any manufactured single component or multiplecomponents. The storage medium or media can be located either in themachine running the machine-readable instructions, or located at aremote site from which machine-readable instructions can be downloadedover a network for execution.

In some embodiments, computing system 1900 contains one or morecompensation module(s) 1908. In the example of computing system 1900,computer system 1901A includes the compensation module 1908. In someembodiments, a single compensation module may be used to perform some orall aspects of one or more embodiments of the methods 400, 500, 1800. Inanother embodiment, a plurality of compensation modules may be used toperform some or all aspects of methods 400, 500, 1800.

It should be appreciated that computing system 1900 is one example of acomputing system, and that computing system 1900 may have more or fewercomponents than shown, may combine additional components not depicted inthe example embodiment of FIG. 19, and/or computing system 1900 may havea different configuration or arrangement of the components depicted inFIG. 19. The various components shown in FIG. 19 may be implemented inhardware, software, or a combination of both hardware and software,including one or more signal processing and/or application specificintegrated circuits.

Further, the steps in the processing methods described herein may beimplemented by running one or more functional modules in informationprocessing apparatus such as general purpose processors or applicationspecific chips, such as ASICs, FPGAs, PLDs, or other appropriatedevices. These modules, combinations of these modules, and/or theircombination with general hardware are all included within the scope ofprotection of the invention.

Geologic interpretations, models and/or other interpretation aids may berefined in an iterative fashion; this concept is applicable to methods400, 500, 1800 as discussed herein. This can include use of feedbackloops executed on an algorithmic basis, such as at a computing device(e.g., computing system 1900, FIG. 19), and/or through manual control bya user who may make determinations regarding whether a given step,action, template, model, or set of curves has become sufficientlyaccurate for the evaluation of the subsurface three-dimensional geologicformation under consideration.

The foregoing description, for purpose of explanation, has beendescribed with reference to specific embodiments. However, theillustrative discussions above are not intended to be exhaustive or tolimit the invention to the precise forms disclosed. Many modificationsand variations are possible in view of the above teachings. Moreover,the order in which the elements of the methods 400, 500, 1800 areillustrate and described may be re-arranged, and/or two or more elementsmay occur simultaneously. The embodiments were chosen and described inorder to best explain the principals of the invention and its practicalapplications, to thereby enable others skilled in the art to bestutilize the invention and various embodiments with various modificationsas are suited to the particular use contemplated.

What is claimed is:
 1. A method for compensating for spatial andslowness or angle blurring of plane-wave reflection coefficients inimaging, comprising: determining one or more source wave fields at areference depth proximate to a reflector for one or more shot records;determining one or more receiver-side blurring functions at thereference depth; constructing an aggregate blurring function based atleast partially on the one or more source wave fields and the one ormore receiver-side blurring functions; and determining, using aprocessor, one or more plane-wave reflection coefficients based at leastpartially on the aggregate blurring function.
 2. The method of claim 1,wherein the one or more plane-wave reflection coefficients are furtherbased at least partially on an interaction of the aggregate blurringfunction with a reflection operator containing the plane-wave reflectioncoefficients.
 3. The method of claim 1, wherein the one or moreplane-wave reflection coefficients are determined based at leastpartially on the aggregate blurring function in a space-time domain at afixed depth proximate to the reference depth.
 4. The method of claim 1,wherein the one or more plane-wave reflection coefficients aredetermined in a domain of lateral position and depth at a fixed time. 5.The method of claim 1, further comprising constructing a calibratedextended image gather proximate to the reference depth based at leastpartially on migration data and a calibration field comprising a bandlimit of the one or more receiver-side blurring functions.
 6. The methodof claim 5, wherein the reflector is a dipping reflector with respect toa coordinate frame of the one or more source wave fields and thecalibrated extended image gather.
 7. The method of claim 6, furthercomprising: transforming the calibrated extended image gather and theaggregate blurring function to a space-frequency domain at a fixeddepth; and obtaining a reflection operator in the space-frequency domainby matrix inversion.
 8. The method of claim 7, further comprisingtransforming the reflection operator from the space-frequency domaininto the one or more plane-wave reflection coefficients by a Fouriertransform.
 9. The method of claim 5, further comprising transforming theaggregate blurring function into a smearing function in a slowness orangle domain of the one or more plane-wave reflection coefficients. 10.The method of claim 9, wherein the smearing function in the slowness orangle domain is determined in a space-time domain at a fixed depthproximate to the reference depth.
 11. The method of claim 9, wherein thesmearing function in the slowness or angle domain is determined in adomain of lateral position and depth at a fixed time.
 12. The method ofclaim 9, further comprising: transforming the calibrated extended imagegather into the slowness or angle domain; and determining the one ormore plane-wave reflection coefficients based at least partially on thesmearing function in the slowness or angle domain.
 13. The method ofclaim 12, further comprising obtaining the calibrated extended imagegather in the smearing function in the slowness or angle domain with aRadon transform or source-direction gathers.
 14. The method of claim 1,further comprising displaying the aggregate blurring function, a rate ofchange in a shape of the aggregate blurring function, the one or moreplane-wave reflection coefficients, or a combination thereof.
 15. Themethod of claim 1, wherein the one or more shot records are generated byone or more shots fired by a user on land or in a marine environment.16. A non-transitory computer-readable medium storing instructions that,when executed by at least one processor of a computing system, cause thecomputing system to perform operations, the operations comprising:determining one or more source wave fields at a reference depthproximate to a reflector for one or more shot records; determining oneor more receiver-side blurring functions at the reference depth;constructing an aggregate blurring function based at least partially onthe one or more source wave fields and the one or more receiver-sideblurring functions; and determining, using the at least one processor,one or more plane-wave reflection coefficients based at least partiallyon the aggregate blurring function.
 17. The non-transitorycomputer-readable medium of claim 16, wherein the one or more plane-wavereflection coefficients are further based at least partially on aninteraction of the aggregate blurring function with a reflectionoperator containing the plane-wave reflection coefficients.
 18. Thenon-transitory computer-readable medium of claim 16, wherein the one ormore plane-wave reflection coefficients are determined based at leastpartially on the aggregate blurring function in a space-time domain at afixed depth proximate to the reference depth.
 19. The non-transitorycomputer-readable medium of claim 16, further comprising constructing acalibrated extended image gather proximate to the reference depth basedat least partially on migration data and a calibration field comprisinga band limit of the one or more receiver-side blurring functions.
 20. Acomputing system, comprising: one or more processors; and a memorysystem comprising one or more non-transitory, computer-readable mediacomprising instructions that, when executed by at least one of the oneor more processors, cause the computing system to perform operations,the operations comprising: determining one or more source wave fields ata reference depth proximate to a reflector for one or more shot records;determining one or more receiver-side blurring functions at thereference depth; constructing an aggregate blurring function based atleast partially on the one or more source wave fields and the one ormore receiver-side blurring functions; and determining, using the one ormore processors, one or more plane-wave reflection coefficients based atleast partially on the aggregate blurring function.